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ABSTRACT 

We apply time- distance helioseismology, local correlation tracking and Fourier 
spatial-temporal filtering methods to realistic supergranule scale simulations of 
solar convection and compare the results with high-resolution observations from 
the SOHO Michelson Doppler Imager (MDI). Our objective is to investigate 
the surface and sub-surface convective structures and test helioseismic measure- 
ments. The size and grid of the computational domain are sufficient to resolve 
various convective scales from granulation to supergranulation. The spatial veloc- 
ity spectrum is approximately a power law for scales larger than granules, with 
a continuous decrease in velocity amplitude with increasing size. Aside from 
granulation no special scales exist, although a small enhancement in power at 
supergranulation scales can be seen. We calculate the time-distance diagram for 
/- and p-modes and show that it is consistent with the SOHO/MDI observations. 
From the simulation data we calculate travel time maps for surface gravity waves 
(/-mode). We also apply correlation tracking to the simulated vertical velocity 
in the photosphere to calculate the corresponding horizontal flows. We compare 
both of these to the actual large-scale (filtered) simulation velocities. All three 
methods reveal similar large scale convective patterns and provide an initial test 
of time-distance methods. 



Subject headings: convection — Sun: oscillations — methods: numerical 
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Introduction 



Local helioseismology uses the observed oscillation signals to probe the near surface 
structure of the Sun to determine sound speed, flow velocities and magnetic field structures. 
One of several methods used in this field is called the time-distance helioseismology. It 
measures th e time taken by a wave packet to travel from one point on the solar surface 
to another (IDuvall et al.l Il993l : iKosovichev fc Duvalll 119971 ). The travel time for a wave 
depends on the wave speed and the flow velocities along the ray paths. These effects can be 
separated by measuring the wave travel time for waves propagating in opposite directions 
along the same ray paths. Inversion methods applied to infer the properties of subsurface 
convection and stru ctures are based on various approxirnation s, such as a geometric acoustic 



ray approxirnation ( Kosovichev et a 



( jjensen et al.l 12001 



Couvidat et al. 



2000l : IZhao et al.l 1200 ll ) , Fresn el-zone approximatio n 



20041 ) and Born-approximation (iGizon &: BirchI 120021). 



Unfo r tunately, until now , there has been no direct verification of these methods (jJensen et al. 
2003l : IWerne et allbooi ). 



Realistic 3D simulations of solar convection ( Stein fc Nordlundll2000l : iBenson et al. hood ) 
provide a way to evaluate the accuracy and consistency of time-distance and other local 
helioseismology techniques. These simulations provide a model of convective motions in the 
upper convection zone and possesses a rich spectrum of /- and p-mode oscillations excited by 
the turbulent convection. They can be used for testing helioseismology methods. Therefore, 
it is very important to apply some of the existing local helioseismology methods to the 
simulated solar convection, and check whether the results are consistent with the subsurface 
properties present in the simulated flows. In addition, the comparison of the oscillation 
properties obtained from the numerical model, such as the velocity spectrum, k — uj diagram 
and time-distance ( cross- covariance) diagram, with observations provides validation for the 
simulation model. 

Previously, the realistic simulations of solar co nvection have be en used to study the 
excitation mechanism of solar and stellar oscillatio ns (IStein et al.ll2004f). an d some properties 



of oscillation modes, such as the line asymmetry (IGeorgobiani et al. 



2000). However, these 



simulations were carried out for small, granulation-size, computational domains. Recently, 
with the progress of parallel supercomputing it became possible to substantially expand the 
computational domain, both horizontally and in dept h, and to simulate the multi-scale solar 
convection, from granulation to supergranular scales (IBenson et al.ll2006l ). These simulations 
allow us for the first time to investigate the properties of large-scale convection by applying 
helioseismology and other methods used to analyze solar observations, and thus test these 
methods and also evaluate how close the realistic simulations are to the real Sun. Such 
tests are particularly important for the local helioseismology methods, which are based on 



- 3 - 



simplified models of wave propagation. 



The goal of this paper is to investigate the basic helioseismic properties of the large-scale 
3D simulations (the oscillation power spectrum and wave travel times) and co mpare these 
with the high- resolution observations from SOHO/MDI (jScherrer et al.lll995l ). An initial 
time-distance analysis is carried out for surface gravity waves (/-mode), which are well suited 
for stu dying the surface and sub-surface structure of solar convection on supergranulation 
scales dPuvall fc Giro^boool ). Our analysis shows that travel times from the simulation 
agree very well with travel times from the SOHO/MDI observations. This demonstrates 
that the simulated data are sufficiently close to the observed solar data, providing ground 
for further detailed helioseismic measurements and inversions. Large-scale structures can be 
detected in the simulated data set using /-mode time-distance and local correlation tracking 
techniques and are observed directly in filtered flow flelds. 



2. Numerical Model and Observed Data 

We study the results of a 3D, compressible, radiative-hydrodynamic (RHD) code simu- 
lating the upper solar convection zone and photosphere. The code calculates LTE, non-gray 
radiation transfer and employs a realistic equat io n of state and opaci ties. More description 



of the code is found in IStein fc Nordlundl (120001 ): iBenson et al.l (120061 ) 



The computational domain of our simulations spans 48 Mm x 48 Mm horizontally and 
20 Mm vertically, with a horizontal resolution of 100 km and vertical resolution of 12 to 75 
km. The three velocity components measured at 200 km above Tcont = 1 are saved every 
minute. Their spatial grid is 500 by 500 pixels. The domain is sufficiently large, and the 
time sequence of several hours is long enough to obtain sufficient signal/noise by temporal 
averaging to perform a helioseismic time-distance analysis and seek evidence for the presence 
of large-scale flows. Eventually, the depth dependence and dynamics of these large-scale 
flows will be studied as well. We compare some of our flndings with results obtained from 
SOHO/MDI high resolution observations. This observational data set is an 8.5 hour time 
series, with 1 minute cadence, of a 211.5 Mm x 211.5 Mm horizontal patch of the MDI 
Doppler velocity, on the spatial grid of 512 by 512 pixels with the resolution of about 400 
km per pixel. The solar rotation is removed. Both sets of data are processed identically, 
whenever a comparison between the observations and simulations occurs. 



-4- 



Surface Structures 



Oscillation modes in the simulations are excited naturally, due to convective motions 



and realistic coo ling at the solar surface (iGeorgobiani et al.l I2OOOI : IStein fc Nordlundl 12001 



Stein et al.l |2004| ). We use the vertical velocity component from the simulation data as a 
proxy for the observed Doppler velocity, and calculate the power spectrum, P{k,uj), in the 
horizontal wavenumber-frequency domain. Following the helioseismology convention, we 
represent the power spectrum in terms of the spherical harmonic degree, £ = kRQ, where 
Rq is the solar radius, and the cyclic frequency, u = ujII-k. 

The power spectrum {Ji — v diagram) calculated from the simulated vertical velocity 
clearly shows a resolved /-mode along with a p-mode spectrum that looks very similar to 
the results from the SOHO/MDI Doppler images (Fig. [1]). These simulations possess a richer 
mode spectrum than in our smaller (6 Mm wide by 3 Mm deep) simulations because the 
wider and deeper domain encompasses many more resonant modes within it. The spec- 
tral resolution achieved in these s i mulat ions allows us to test time-distance helioseismology 
measurements (iDuvall et al.lll993l . 119971 ). 



Mode filtering is an essential part of the time-distance technique. The convection sig- 
nals (broad wedge in the lower right part of the £ — v diagram) need to be removed. We 
construct the time-distance diagram for both simulated and observed data b y computing 



cross -covariance of oscillation signals separated by a certain distance (Fig. [2]) (IDuvall et al. 



19971 ). and find a remarkable similarity between them. Because of the higher resolution, 
in the numerical simulation the cross-covariance signal extends to much shorter distances 
than in the MDI data. This illustrates a potential for small-scale helioseismic diagnostics 
of observati ons of higher than MDI resolution, such as anticipated from the Solar-B mission 
JSekiill2004l ). 



We show in section 13.21 that the travel times from simulation data agree with travel 
times from the SOHO/MDI observations very well. Both the i — v and time-distance di- 
agrams illustrate the potential of the simulations for testing helioseismology measurement 
and inversions procedures. 

In the following sections, we explore the surface structures in the simulated data, directly 
by averaging the flow field and using spectral analysis (Sect. 13.11) . and also by implement- 
ing the /-mode time-distance analysis (Sect. 13.21) . and local correlation tracking technique 
(Sect.[33D. 
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Spherical harmonic degree Spherical harmonic degree 



Fig. 1. — The power spectra {i — u diagrams) for the simulated vertical velocity (left) and 
the Doppler velocity from the MDI high- resolution observations (right). The dark curve 
represents the theoretical /-mode. 



3.1. Fourier Analysis of Surface Velocities 



We analyze the simulated velocity fields by calculating spatial power spectra for both 
vertical and horizontal velocity components and compare these with the spatial power spec- 
trum for the MDI Doppler signal. Power spectra are calculated for each time and are then 
averaged over time. This reduces the statistical fluctuations that are present in spectra from 
a single time. In Fig. [3] we show V{i) = (£P(^))^/^, where V{i) is velocity amplitude in the 
Fourier domain at spherical harmonic degree £ and P{i) is time averaged power (velocity 
squared) per unit i. There is a continuous decrease in the velocity amplitude from granules 
(~ 1.5 Mm, i ~ 2000-3000) to larger scales. Aside from granulation no special scales are re- 
vealed in the spectra, although there may be a small enhancement in the horizontal velocity 
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Fig. 2. — Time- distance diagrams for the simulated vertical velocity (left) and the Doppler 
velocity from the high-resolution MDI observations (right). The dark curve is time-distance 
relation computed from a standard solar model. The second bounce can be seen in both 
pictures. 

(^oriz) spectrum at supergranulation scales (£ ~ 200). 

In Fig. m we separate motions into 'convective' (dashed lines) and 'oscillatory' (dotted 
lines) parts, according to whether the power in Fourier space lies above or below the line 
uj = ck, where c = 6 kms~^ is approximately equal to the sound speed at the photosphere. 
As can be seen from Fig. [T] this provides a rather clean split; the corresponding line in Fig. [T] 
would intersect the right hand side u-axis at z/ 3.8 mHz. 

As shown by Fig. |H the horizontal velocities (blue curves) at the height of the MDI 
Ni A 676.78 nm line formation (about 200 km above t^qq = 1) are almost exclusively of 
convective origin. The vertical velocities (green curves), on the other hand, are a mixture of 
convective (dashed line) and oscillation (dotted line) motions where the convective motions 
are dominant at smaller scales and the oscillatory motions are dominant at larger scales. 

For granular and larger scales the horizontal velocities are larger than the vertical ve- 
locities and become increasingly dominant as the scale increases. The reason that large 
scale horizontal motions are present at these heights is that the atmosphere is not far from 




Fig. 3. — Spatial spectra for horizontal (dotted) and vertical (dashed) components of the 
simulated velocity, and the MDI Doppler signal (solid). Spectra are calculated for individual 
snapshots and then time-averaged, the simulation data over 20 hours and the MDI data over 
8.5 hours. 

hydrostatic equilibrium and the pressure fluctuations that exist below the surface to drive 
the large scale horizontal flows imprint their fluct uations on the s urface and produce smaller 



amplitude large scale horizontal flows there also (jNordlundlll982l ). 



Because of mass conservation, vertical convective velocities decrease more rapidly with 
size than the horizontal ones, with the ratio of vertical to horizontal velocity amplitudes 
approximately inversely proportional to size (cf. the green and blue dashed lines in Fig. H]). 
At granular scales, where the vertical velocity peaks, the vertical convective velocity at the 
height of formation of the A 676.78 nm line is still about 3-4 times weaker than the horizontal 
one. 




10^ 10^ 10^ 10^ 10^ 

Spherical harmonic degree 



Fig. 4. — Velocity spectra for the simulation (vertical = green, horizontal = blue) and MDI 
observations (red) separated into convective (dashed lines), oscillatory (dotted lines) and 
total (solid lines) components. The oscillatory signal is negligible in the horizontal velocity 
but dominates the vertical velocity and MDI Doppler velocity on large scales < 1000). 

In general, the MDI Doppler signal is a mix of horizontal and vertical velocities. There is 
no obvious way to decouple these components, although, depending on where on the solar disc 
the data patch was located, one can estimate the approximate contribution of the horizontal 
and vertical components. The current dataset was taken near solar disk center (which at the 
time was about 7.5 heliographic degrees below the solar equator). Nevertheless, because of 
the significant extension of the patch, the root-mean-square projection of horizontal velocities 
into the line-of-sight was about 15%, while almost 99% of the vertical velocity is projected 
onto the line-of-sight. Indeed, if we combine 15% of the simulated horizontal velocities and 
99% of the simul ated vertical velocitie s, multiply by the instrument Modulation Transfer 



Function (MTF) (IWoodard et al.ll200ll ) and compare the resulting power spectrum with the 
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MDI power spectrum in Fig. [5], we get a good correspondence, although the simulated power 
at large scales (mostly oscillatory according to Fig. H]) falls somewhat short of the observed 
power. Given that the numerical res olution (96 km) is somew hat marginal on granular scales, 
where solar oscillations are excited (IStein fc Nordlundll200ll ). this is not surprising. 



From the slope of the MDI convective velocity (red dashed) curve in Fig. H] one surmises 
that the convective LOS velocity observed with MDI over this patch actually derives most 
of its contribution at large scales from the horizontal velocity field. An enhancement at 
supergranular scales (15-55 Mm, ^ ~ 80-300) is visible. 




Fig. 5. — Simulated line-of-sight (LOS) component of the surface velocity spectrum (dotted), 
convolved with the instrumental MTF (dashed), compared with the MDI high- resolution 
spectrum (solid). 

Averaging the velocity fields over several hours eliminates the oscillations and suppresses 
the small scale features but does not eliminate them. To remove the small spatial scales. 
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it is necessary, in addition, to apply a low-pass k-space filter, where all values outside the 
circle < 0.1 Mm~^ are set to zero. We use this procedure to reveal the large- 

scale patterns in both vertical and horizontal velocity components. Fig. [6] (left panel) is an 
example of such a pattern in the component. An alternative, cleaner, procedure would be 
to filter entirely in Fourier space, eliminating the oscillations by removing power for u > ck 
(as discussed in connection with Fig. H] above) and then applying a low-pass filter in k. Here 
we follow the first procedure which has been generally used in time-distance helioseismology. 



Vz Power map 




X (Mm) X (Mm) 



Fig. 6. — Left: Simulated large scale vertical velocity. The velocity is averaged over 8.5 
hours and smoothed with a low-pass filter. Light is downflow and dark upflow. Right: p- 
mode power in the frequency interval 2-5 mHz. No low-pass filter has been applied. The 
oscillatory power is concentrated in the intergranular lanes and is weakest in the upflow 
centers of the large-scale structures. 

We also construct p-mode power maps by Fourier-transforming the vertical velocity 
fields in time at each spatial point, and summing the result over the frequency interval 2-5 
mHz (Fig. El right panel). The p-mode power is rather i ntermittent on sm a ll sca les, and is 



in fact concentrated in intergranular lanes, as shown by IStein fc Nordlundl (120011 ). A large 
scale modulation, correlated with the large scale vertical velocity field, is also discernible. 
Presumably the large scale modulation is caused by the slightly larger velocity amplitudes 
present in intergranular lanes inside supergranulation scale downdrafts. 
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3.2. F-Mode Time-Distance Analysis 

An important question is wlietlier local helioseismology methods can recover the flow 
pattern from the oscillatory component of the surface velocity fleld in these simulations. 
We applied the time-distance technique to the /-mode oscillation, because /-mode is the 
surface gravity mode and thus contains information about horizontal flows and structures 
near the surface. Therefore, the time-distance results can be directly compared with the 
surface properties. The /-mode analysis is commonly used in helioseismology. We isolate the 
/-mode in kx, ky, u space by retaining only values along a band centered on the theoretical /- 
mode ridge and setting all other values to zero, and then perform the time-distance analysis, 
constructing outgoing and ingoing travel time maps, mean travel times and travel time 
differences. We calculate the horizontal divergence dVx/dx + dVy/dy from the simulation 
data and co mpare it to the travel t ime differences map after a low-pass flltering (Fig. [7]). 



As shown by iDuvall fc GizonI (120001 ). the travel-time difference is roughly proportional to 



the horizontal flow divergence. We see a good agreement in the large-scale structures, with 



- divergence oi travel time 




12 24 36 48 12 24 36 48 

X (Mm) X (Mm) 



Fig. 7. — Comparison of the simulated horizontal velocity divergence (left; converging flows 
are light, diverging flows are dark) with /-mode outgoing and ingoing travel time difference 
map (right; incoming waves are faster in light areas, while outgoing waves are faster in dark 
areas) after a low-pass fllter (as in Fig. [6]) has been applied. Travel time extrema are 50 s 
for outgoing and 73 sec for incoming waves. 
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correlation coefficient 0.88. Calculating east-west and north-south travel time differences, 
we obtain proxies for horizontal velocity components Vx and Vy. We compare these time 
difference maps with the actual velocities (Fig. [H] for the North-South componenet, the 
other looks similar). After low-pass filtering, there is a good qualitative agreement between 
the two, with correlation coefficients 0.7 for x-component and 0.73 for y-component. It is 



VY 



-ns 
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Fig. 8. — Comparison between the simulated horizontal velocity {Vy, left, light is north- 
bound and dark is south-bound) and /-mode north-south travel time difference map (right, 
light represents faster north-bound and dark represents faster south-bound waves). A low- 
pass filtering has been performed here as well. Travel time extrema are 39 s for north-bound 
and 72 s for south-bound waves. 

worth mentioning that the time-distance technique measures an average over a certain depth, 
therefore a very high correlation is not to be expected when comparing the time-distance 
results to a single depth data, particularly above the surface. The agreement might have 
been better if one looked at a depth representative of the /-mode travel paths (/-mode 
kinetic energy is concentrated in the 2 Mm just below the surface). 
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3.3. Local Correlation Tracking 



Another possibility for deducing the horizontal velocities is to apply the local corre- 
lation tracki ng (LCT) technique to th e vertical velocity. The LCT is a cross-correlation 



method (e.g. [November fc Simonlll988l ) applied to a time series of solar granulation images. 



The cross correlation is spatially localized (within a certain window, usually a Gaussian); 
its time average presumably measures horizontal displacements of the flows. We apply the 
LCT technique to the simulated vertical velocity field and compare the resulting horizontal 
velocity proxies with the simulated horizontal velocity components. Fig. [9] shows the Vy 
component. Velocities determined by the LCT method are tightly correlated with the simu- 
lation velocity, with correlation coefficients of 0.99 for both and Vy components. However, 
the LCT amplitudes are a factor 1.8 smaller. Application of a low-pass filter to the simu- 
la tion velocities b r ings t his ratio down to 1.5. These results are consistent with conclusions 
of iRieutord et al.l (120011 ). According to them, the horizontal velocity proxies deduced from 
granular motions underestimate actual velocities by a variable factor, from 2.1 at small scales 
to 1.6 at large scales. 



VY (SIM) 
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Fig. 9. — Simulated horizontal velocity component Vy (left), averaged over 8.5 hours, and 
Vy, obtained by LCT analysis (right), after low-pass filtering. 

Thus, the realistic simulation data also provide a tool for testing and improving the 
local correlation tracking techniques. 
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4. Discussion 

The new large-scale realistic 3D radiative-hydrodynamic simulations of the upper layers 
of solar convection provide an excellent opportunity to validate various techniques, widely 
used in solar physics and helioseismology for directly obtaining otherwise inaccessible prop- 
erties (subsurface flows and structures, etc). On the other hand, these analysis techniques 
also help to examine how realistic the simulations are. We have performed an initial anal- 
ysis of the simulated data and compared our results with the outcome of the SOHO/MDI 
observations. The similarity between the simulated and observed I — v and time-distance 
diagrams demonstrates that the simulations can be efficiently used to perform and validate 
local helioseismology techniques. This agreement also reveals a potential for high-resolution 
time-distance measurements. We carried out /-mode time-distance calculations and local 
correlation tracking to obtain horizontal velocities. We compare them to the simulated hor- 
izontal flows and flnd a good qualitative agreement. We see large-scale structures in both 
actual horizontal velocities and their proxies. The results of this investigation provide the ba- 
sis for further detailed helioseismic analyses of the simulated data, including various scheme 
of measuring p-mode travel times, approximations for the sensitivity kernels and inversion 
procedures. 
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